Fast realization of a spatially correlated percolation model 



School of Science, 



Hongting Yang* 
Wuhan University of Technology, 



Wuhan 430070, P.R. China 



o 

(N 
C 

o 



i 

C/3 



c3 



i 

C 

O 

o 



> 
o 

o 



Stephan Haas 

Department of Physics and Astronomy, University of Southern California, Los Angeles, CA 90089-0484 

(Dated: January 31, 2013) 

We propose two schemes to achieve fast realizations of spatially correlated percolation models. 
The schemes are shown to be efficient in complementary regimes of correlation phase space. They 
are combined with a generalized Newman-Ziff algorithm to numerically determine the percolation 
thresholds of two-dimensional lattices in the presence of correlations, ft is found that the spatial 
correlations affect only a relatively small part of phase space. 
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I. INTRODUCTION 

Much of the existing literature on percolation has fo- 
cused on the uncorrelated case where each lattice 
site or bond is independently occupied with probability 
p or empty with probability 1 — p. Some of the occu- 
pied sites/bonds form clusters, and with more sites or 
bonds occupied, these clusters grow larger. For finite lat- 
tices with periodic boundary conditions, if a cluster grows 
large enough to wrap around the lattice, wrapping perco- 
lation occurs [2J . Analogously, for finite lattices with free 
boundaries, the appearance of a spanning cluster touch- 
ing two opposite boundaries indicates the occurrence of 
spanning percolation along the direction perpendicular 
to the two boundaries. 

The correspondence of phase transition with perco- 
lation has attracted much interest [3, Ej|. Applications 
range widely from quarks and gluons to galaxies [H-fl], 
and from epidemics to networks [9l4l4|. In particular, 
the notion of discontinuity in explosive percolation has 
stimulated a number of interesting recent studies 
[l9j . For the square lattice, it is believed that the value 
of the percolation threshold is unique, although its exact 
value remains unknown for site percolation [20l422l |. 

While most of these previous studies do not include 
spatial correlations, this constraint appears overly re- 
strictive. For example, phase transitions between solid, 
liquid and gaseous states take place under different con- 
ditions. A theory with unique percolation threshold can 
not be applied to describe multiple phase transitions in 
a unified way. Effective percolation models without in- 
teractions between sites or bonds are incomplete [l2| . 
Such interactions, whether attractive or repulsive, gener- 
ate spatially correlated site and bond distributions. How- 
ever, in contrast to uncorrelated percolation models, it is 
more complicated to efficiently generate spatially corre- 
lated systems [23j. 

In this brief report, we examine percolation models 
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with spatial correlations caused by attractions between 
occupied sites. We will address two issues. The first 
question is how to generate spatially correlated distribu- 
tions efficiently. In the following section, we present two 
schemes which work well in complementary parameter 
regimes. The second question is how to find the value of 
a percolation threshold. A generalized Newman-Ziff al- 
gorithm provides the tool for this second task 24] . This 
will be discussed in the third section. 



II. MODEL AND ALGORITHM 

Here we discuss how to build spatially correlated lat- 
tices by considering compactness among occupied sites. 
It is assumed that there is a bond between all neigh- 
boring sites. In terms of the bond length (the lattice 
spacing), the least number of bonds traversed to reach 
one site from another site is defined as the distance be- 
tween the two sites. If two sites satisfy the correlation 
condition, i.e. the distance between the two sites is no 
greater than a preset value, we say that the two sites 
are spatially correlated. We use d to represent this pre- 
set value, and hence the correlation condition is given 
by d = 1,2,3, ... Obviously, d — 1 corresponds to 
the most compact distribution, and d = L represents a 
completely uncorrelated random distribution. 

Let us now contemplate algorithms how to generate 
spatially correlated distributions for a square lattice with 
linear dimension L. Only the first occupied site is chosen 
completely randomly. The distance between the next site 
to be occupied and at least one of the already occupied 
sites is set to be no greater than d. This way, each site 
occupied later is spatially correlated to (at least) one of 
the previously occupied sites. For d = 1, since this is 
the most compact distribution, the effective attraction 
between the occupied sites is the strongest, and in the 
process of occupying sites one by one, there is only one 
cluster. Two examples of spatially correlated distribu- 
tions for d — 1 and d — 2 on square lattices with L = 8 
are shown in Fig. 1(a) and l(b)| 

We now attempt to answer the question of how to gen- 
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(a) d = 1 



(b) 



FIG. 1: Two examples of spatially correlated distributions 
on a periodic square lattice with L = 8. Ten occupied sites 
shaded in (a) and (b) are shown for d = 1 and d — 2 respec- 
tively. 



erate a spatially correlated percolation model most effi- 
ciently. While this may not be a pressing issue for small 
lattices as those shown in Fig. [TJ for large lattices with L 
up to 256, inefficient methods to generate spatially corre- 
lated distributions of occupied sites will take exceedingly 
long running time. The key to numerical efficiency is 
when we should test the correlation condition. In gen- 
eral, there are two schemes, which are discussed below. 

In scheme A, we test the correlation condition after 
we randomly choose an empty site, which is chosen to be 
occupied if it satisfies the correlation condition. Other- 
wise, we continue to randomly choose other empty sites, 
until we find one that fulfills the correlation condition. 
Obviously, for a randomly chosen site, greater d implies 
a greater probability to satisfy the correlation condition. 
In contrast, a smaller d means that the empty sites al- 
lowed by the correlation condition are distributed in a 
smaller area around the occupied sites, and it will there- 
fore take a longer time to find a correlated empty site 
among all the empty sites. Scheme A can be summa- 
rized as follows. 

(1) Randomly choose one empty site. 

(2) If this site is spatially correlated to one of the pre- 

viously occupied sites, occupy this site; otherwise, 
go to step (1). 

Repeat steps (1) to (2) until all sites are occupied or 
wrapping percolation occurs. 

In scheme B, the correlation condition is tested before 
we randomly choose an empty site. We test the correla- 
tion condition, and find all available empty sites spatially 
correlated to the newly occupied site. We then randomly 
choose one empty site to occupy directly from this list of 
available empty sites. This scheme works effectively for 
the case of small d on large lattices, whereas for large d on 
a large lattice, it will take a longer time to test the corre- 
lation condition and find the spatially correlated empty 
sites. Scheme B can be summarized as follows. 

(1) Occupy one site randomly chosen from the list of cor- 
related empty sites, which are spatially correlated 
to the previously occupied sites. 



(2) Refresh the list of correlated empty sites by adding 
extra empty sites spatially correlated to the newly 
occupied site in step (1). 

These steps are repeated until all sites are occupied or 
wrapping percolation occurs. 

To show the great difference between the efficiency of 
the two schemes in calculation time, an example of a 
lattice L = 256 is considered in Table Q] The number of 
runs of the algorithm for each scheme is taken as 2 x 10 4 . 
The computations are implemented on a desktop PC with 
a CPU clock speed 2.6 GHz and memory 1.96 GB. When 
d = 1, scheme B takes about 336 seconds, while scheme 
A takes about 260 days. When d < 16, scheme B takes 
shorter time than scheme A. When d > 32, the former- 
takes longer time than the latter; when d — 128, the 
former takes about 65 hours while the latter takes only 
260 seconds. Obviously, scheme A is not effective when 
d is small, and scheme B works less effectively when d is 
large. For the case of a completely uncorrelated random 
distribution, d = L, scheme B will take about 8 days, 
while scheme A takes about 210 seconds, which is close to 
170 seconds spent in the Newman-Ziff algorithm. Similar 
results exist for other lattices with different L. With 
increasing L, the difference of the two schemes in the 
calculation time increases too. 

In general, in scheme B, the computation time Tl ~ 
nL 2 for fixed d/L > |; for fixed L, Tl increases with in- 
creasing d, Tl ~ d Cl , where c\ — 1.1, 1.3, 1.4, 1.4 for L = 
32, 64, 128, 256 respectively. In scheme A, the computa- 
tion time Tl ~ nL 2 for fixed d; for fixed L, Tl decreases 
with increasing d, Tl ~ d~ c ' 2 , where C2 = 1.7, 2.0, 2.1, 2.3 
for L = 32, 64, 128, 256 respectively. Obviously, from the 
point of saving computation time, any scheme alone is 
not appropriate to the calculation for the whole range 
of the values of d's. The two schemes should be com- 
bined to give a fast realization of the spatially correlated 
percolation model. 



III. EXTRACTION OF PERCOLATION 
THRESHOLD 

On a lattice with N sites, the Newman-Ziff algorithm 
states that a state with n + 1 occupied sites is achieved 
by adding one extra randomly chosen site to a state with 
n occupied sites. An observable Q(p) is binomially ex- 
panded in terms of the measured quantities {Q n } 



(1) 



Such observables Q can be the probability of cluster- 
wrapping, mean cluster size, correlation length, and so 
on. Spatial correlations affect the values of Q n , and thus 
the value of Q(p). There are four types of probabilities 
Rl (p) of cluster wrapping on the periodic square lattice 
of A = LxL sites. R L is the probability of cluster wrap- 
ping along either the horizontal or vertical directions, or 
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TABLE I: Comparison of calculation time for the two schemes at different d's, on the lattice with L — 256. The number of 
runs of the algorithm for each scheme is n — 2 x 10 4 . 



d 


1 


2 


4 


8 


16 


32 


64 


128 


256 


B-scheme 


336 s 


492 s 


885 s 


30 min 


85 min 


4.7 h 


17 h 


65 h 


8 d 


W-scheme 


260 d 


60 d 


11 d 


42 h 


6.4 h 


1 h 


700 s 


260 s 


210 s 



Here, d=day(s), h=hour(s), min=minute(s), s=second(s). 



both; R L , around one specified axis but not the other 
axis; , around both the horizontal and vertical direc- 
tions; and R£ , around the horizontal and vertical 
directions, respectively. The four wrapping probabilities 
satisfy the equations 
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only two of which need to be measured. In general, given 
the exact value of Roo(p c ), the solution p of the equation 



Rl(p) = RooiPc), 



(4) 



gives a very good estimator for the threshold p c in per- 
colation theory. However, for spatially correlated perco- 
lation models, the exact value of Roo(Pc) is unknown in 
advance, so we can not obtain the value of p c by solving 
the equation. However, we can obtain p c from Rv' (p) 
by virtue of its non-monotonicity. We use Machta's 
method 0, [H| , instead of the alternative criterion for 
two-dimensional wrapping percolation (22j , to save com- 
putation time. 



than the nearby p c values. With increasing L, there are 
no more signs of dip. Therefore, the dip in the curve for 
L = 32 is believed to be the result of finite size effects. 
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FIG. 2: Variations of the percolation threshold with the cor- 
relation strength for the lattices with L — 32, 64, 128, and 
256. 



IV. RESULTS 

All computations are implemented on the same desk- 
top PC. We take d = 1,2,4,... , L/2, L, and restrict our- 
selves to calculating p c . We choose scheme B when d < 4 
for L = 32, d < 8 for L = 64 and L = 128, d < 16 for 
L = 256. In the other cases, we choose scheme A instead. 
Each time of running the program, the algorithm is im- 
plemented for n runs to output the value of p c , and we 
run the program g times to estimate its error. 

For L = 32, d = 1, we take n = 2 x 10 8 , running the 
program once takes about 13 hours. We repeat this ten 
times, and obtain the value of percolation threshold and 
its error: p c — 0.6026982(32). For all other values of L 
and d, we also take g = 10, and the number of runs of 
the algorithm n ranges from 1.4 x 10 5 to 2 x 10 8 , and 
each run takes 8-13 hours. 

To elucidate the effect of the size of a lattice, the per- 
colation thresholds are shown in Fig. |]as a function 
of (d/L) 2 , instead of d. The errors in p c range from 
8.5 x 10~ 7 to 4.6 x 10~ 5 , which are too small to be 
shown. One notices that there is a sharp dip around 
d = 2 only in the curve for L = 32. In the curve for 
L = 64, p c = 0.5926966(31) at d = 8 is a little bit smaller 



0.592748 
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FIG. 3: The values of percolation thresholds under the ex- 
treme conditions, (a) d = L no spatial correlation, and (b) 
d — 1 the most compact distribution. 

The most striking feature in Fig. [2] is that the effect 
of spatial correlation on p c exists only in a very small 
regime of correlation strength (d/L) 2 < 0.004, estimated 
from the existing data. Except the dip in the curve for 
L = 32, the values of p c nearly keep unchanged in a wide 
range of correlation strength 0.004 < (d/L) 2 < 1. 
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Since d = L corresponds to the case of completely un- 
correlated random distribution, the value of percolation 
threshold converges to p c (oo) according to [24( 

Pc -p c (oo) ~L- n /\ (5) 

As in Fig. 3(a), the finite size scaling of estimates for 
p c gives Pc(oo) = 0.59274640(52), which coincides with 
the previous results [20M22I [24| . The strongest correlated 
case is for d = 1, i.e. the strongest effective attraction be- 
tween occupied sites on an infinite lattice. Linear fitting 
in Fig. 3(b) gives p c (d = 1, L = 00) = 0.701(14). 

V. CONCLUSION 

Scheme A and B are combined to give an efficient re- 
alization of spatially correlated lattices, introduced phc- 
nomenologically by restricting the distance between sites. 



Scheme B is effective especially when there are only a few 
sites can be chosen to occupy in the process of cluster 
growing. For any population of sites on a lattice, one 
has to choose an appropriate algorithm to achieve the 
corresponding population. 

As for a spatially correlated percolation model, the cal- 
culation of critical exponents is of great interest, too. 
For the purposes of showing the differences of the two 
schemes in realizing a site population and the effects of 
correlation strength, here we focused on the calculation 
of percolation thresholds. 

Very strong attractions between occupied sites may se- 
riously hinder the occurrence of percolation transition, 
while correlations which are not very strong between oc- 
cupied sites have less effects on the percolation thresh- 
olds. This is the reason why the uncorrelated percola- 
tion theory without considering spatial correlations has 
so many successful applications. 
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